Lecture

Introduction to Distance sampling with inlabru

Jafet Belmont

Jafet.BelmontOsuna@glasgow.ac.uk

School of Mathematics and Statistics, University of Glasgow

November 13, 2025

Outline

What are INLA and inlabru?
Why the Bayesian framework?
Which model are inlabru-friendly?
What are Latent Gaussian Models?
How are they implemented in inlabru?

What is INLA? What is `inlabru`?

The short answer:

INLA is a fast method to do Bayesian inference with latent Gaussian models and inlabru is an R-package that implements this method with a flexible and simple interface.

The (much) longer answer:

Rue, H., Martino, S. and Chopin, N. (2009), Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71: 319-392.
Van Niekerk, J., Krainski, E., Rustand, D., & Rue, H. (2023). A new avenue for Bayesian inference with INLA. Computational Statistics & Data Analysis, 181, 107692.
Lindgren, F., Bachl, F., Illian, J., Suen, M. H., Rue, H., & Seaton, A. E. (2024). inlabru: software for fitting latent Gaussian models with non-linear predictors. arXiv preprint arXiv:2407.00791.
Lindgren, F., Bolin, D., & Rue, H. (2022). The SPDE approach for Gaussian and non-Gaussian fields: 10 years and still running. Spatial Statistics, 50, 100599.

Where?

Website-tutorials

inlabru https://inlabru-org.github.io/inlabru/

R-INLA https://www.r-inla.org/home

Discussion forums

inlabru https://github.com/inlabru-org/inlabru/discussions

R-INLA https://groups.google.com/g/r-inla-discussion-group

Books

Blangiardo, M., & Cameletti, M. (2015). Spatial and spatio-temporal Bayesian models with R-INLA. John Wiley & Sons.
Gómez-Rubio, V. (2020). Bayesian inference with INLA. Chapman and Hall/CRC.
Krainski, E., Gómez-Rubio, V., Bakka, H., Lenzi, A., Castro-Camilo, D., Simpson, D., … & Rue, H. (2018). Advanced spatial modeling with stochastic partial differential equations using R and INLA. Chapman and Hall/CRC.
Wang, X., Yue, Y. R., & Faraway, J. J. (2018). Bayesian regression modeling with INLA. Chapman and Hall/CRC.

So… Why should you use `inlabru`?

What type of problems can we solve?
What type of models can we use?
When can we use it?

Distance Sampling

Overview of Distance Sampling

Distance sampling is a family of related methods for estimating the abundance and spatial distribution of wild populations.
Distance sampling is based on the idea that animals further away from observers are harder to detect than animals that are nearer.

This idea is implemented in the model as a detection function that depends on distance.
- Species at greater distances are harder to detect and the detection function therefore declines as distance increases.

Density surface models

Coupling distance sampling data with spatial modelling allows maps of spatially varying density to be produced.

Traditionally, this is achieved in a two-stages approach by (i) using a detectability point estimates to create an offset vector to (ii) use within GLM or GAM for count response data.

Density surface models

Coupling distance sampling data with spatial modelling allows maps of spatially varying density to be produced.

Traditionally, this is achieved in a two-stages approach by (i) using a detectability point estimates to create an offset vector to (ii) use within GLM or GAM for count response data.
This requires binning the data into counts based on some discretisation of space.

Density surface models

Coupling distance sampling data with spatial modelling allows maps of spatially varying density to be produced.

Traditionally, this is achieved in a two-stages approach by (i) using a detectability point estimates to create an offset vector to (ii) use within GLM or GAM for count response data.
A major downside to this approach is the propagation of uncertainty from the detection model to the second-stage spatial model.

The goal: one-stage distance sampling model, simultaneously estimating the detectability and the spatial distribution of animals using a point process framework.

Point Processes

Point process data

Many of the ecological and environmental processes of interest can be represented by a spatial point process or can be view as an aggregation of one.

Many contemporary data sources collect georeferenced information about the location where an event has occur (e.g., species occurrence, wildfire, flood events).
This point-based information provides valuable insights into ecosystem dynamics.

Defining a Point Process

Consider a fixed geographical region $A$.
The set of locations at which events occur are denoted by $\mathbf{s} = (\mathbf{s}_1, \ldots, \mathbf{s}_n)$.
We let $N(A)$ be the random variable which represents the total number of events in region $A$.
Our primary interest is in measuring where events occur, so the locations are our data.

Homogeneous Poisson Process

The simplest version of a point process model is the homogeneous Poisson process (HPP).
The likelihood of a point pattern $\mathbf{y} = \left[ \mathbf{s}_1, \ldots, \mathbf{s}_n \right]^\intercal$ distributed as a HPP with intensity $\lambda$ and observation window $\Omega$ is

\[ p(\mathbf{y} | \lambda) \propto \lambda^n e^{ \left( - |\Omega| \lambda \right)} , \]
- $|\Omega|$ is the size of the observation window.
- $\lambda$ is the expected number of points per unit area.
- $|\Omega|\lambda$ the total expected number of points in the observation window.

A key property of a Poisson process is that the number of points within a region $A$ is Poisson distributed with constant rate $|A|\lambda$.

Inhomogeneous Poisson process

The inhomogeneous Poisson process has a spatially varying intensity $\lambda(\mathbf{s})$.

The likelihood in this case is

\[ p(\mathbf{y} | \lambda) \propto \exp \left( -\int_\Omega \lambda(\mathbf{s}) \mathrm{d}\mathbf{s} \right) \prod_{i=1}^n \lambda(\mathbf{s}_i). \]

If the case of an HPP the integral in the likelihood can easily be computed as $\int_\Omega \lambda(\mathbf{s}) \mathrm{d}\mathbf{s} =|\Omega|\lambda$
For IPP, integral in the likelihood has to be approximated numerically as a weighted sum.

Inhomogeneous Poisson process

The integral in is approximated as $\sum_{j=1}^J w_j \lambda(\mathbf{s}_j)$

$w_j$ are the integration weights
$\mathbf{s}_j$ are the quadrature locations.

This serves two purposes:

Approximates the integral
re-write the inhomogeneous Poisson process likelihood as a regular Poisson likelihood.

Inhomogeneous Poisson process

The idea behind the trick to rewrite the approximate likelihood is to introduce a dummy vector $\mathbf{z}$ and an integration weights vector $\mathbf{w}$ of length $J + n$

\[\mathbf{z} = \left[\underbrace{0_1, \ldots,0_J}_\text{quadrature locations}, \underbrace{1_1, \ldots ,1_n}_{\text{data points}} \right]^\intercal\]

\[\mathbf{w} = \left[ \underbrace{w_1, \ldots, w_J}_\text{quadrature locations}, \underbrace{0_1, \ldots, 0_n}_\text{data points} \right]^\intercal\]

Then the approximate likelihood can be written as

\[ \begin{aligned} p(\mathbf{z} | \lambda) &\propto \prod_{i=1}^{J + n} \eta_i^{z_i} \exp\left(-w_i \eta_i \right) \\ \eta_i &= \log\lambda(\mathbf{s}_i) = \mathbf{x}(s)'\beta \end{aligned} \]

This is similar to a product of Poisson distributions with means $\eta_i$, exposures $w_i$ and observations $z_i$.
This is the basis for the implementation of Cox process models in inlabru, which can be specified using family = "cp".

Limitations with IPP

IPP models assume that data points are conditionally independent given the covariates, meaning that any spatial variation is fully explained by environmental and sampling factors.
Unmeasured endogenous and exogenous factors can create spatial
Ignoring them can lead to bias in our conclusions.

The Log-Gaussian Cox Process

Log-Gaussian Cox processes (LGCP) extends the IPP by allowing the intensity function to vary spatially according to a structured spatial random effect.

\[ \log~\lambda(s)= \mathbf{x}(s)'\beta + \xi(s) \]

The events are then assumed to be independent given $\xi(s)$ - a GMRF with Matérn covariance.
inlabru has implemented some integration schemes that are especially well suited to integrating the intensity in models with an SPDE effect.

Gaussian Random Fields

If we have a process that is occurring everywhere in space, it is natural to try to model it using some sort of function.

If $z$ is a vector of observations of $z(\mathbf{s})$ at different locations, we want this to be normally distributed:

\[ \mathbf{z} = (z(\mathbf{s}_1),\ldots,z(\mathbf{s}_m)) \sim \mathcal{N}(0,\Sigma) \]

where $\Sigma_{ij} = \mathrm{Cov}(z(\mathbf{s}_i),z(\mathbf{s}_j))$ is a dense $m \times m$ matrix.

Gaussian Random Fields

A Gaussian random field (GRF) is a collection of random variables where observations occur in a continuous domain, and where every finite collection of random variables has a multivariate normal distribution

Stationary random fields

A GRF is stationary if:

has mean zero.
the covariance between two points depends only on the distance and direction between those points.

It is isotropic if the covariance only depends on the distance between the points.

The SPDE approach

The goal: approximate the GRF using a triangulated mesh via the so-called SPDE approach.

The SPDE approach represents the continuous spatial process as a discretely indexed Gaussian Markov Random Field (GMRF)

We construct an appropriate lower-resolution approximation of the surface by sampling it in a set of well designed points and constructing a piecewise linear interpolant.

The SPDE approach

A GF with Matérn covariance $c_{\nu}(d;\sigma,\rho)$ is a solution to a particular PDE.

\[ c_{\nu}(d;\sigma,\rho) = \sigma^2\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{8\nu}\frac{d}{\rho}\right)^{\nu}K_{\nu}\left(\sqrt{8\nu}\frac{d}{\rho}\right) \]

This solution is then approximated using a finite combination of piecewise linear basis functions defined on a triangulation .
The solution is completely defined by a Gaussian vector of weights (defined on the triangulation vertices) with zero mean and a sparse precision matrix.
How do we choose sensible priors for $\sigma,\rho$?

Penalized Complexity (PC) priors

Penalized Complexity (PC) priors proposed by Simpson et al. (2017) allow us to control the amount of spatial smoothing and avoid overfitting.

PC priors shrink the model towards a simpler baseline unless the data provide strong evidence for a more complex structure.
To define the prior for the marginal precision $\sigma^{-2}$ and the range parameter $\rho$, we use the probability statements:
- Define the prior for the range $\text{Prob}(\rho<\rho_0) = p_{\rho}$
- Define the prior for the range $\text{Prob}(\sigma>\sigma_0) = p_{\sigma}$

Learning about the SPDE approach

F. Lindgren, H. Rue, and J. Lindström. An explicit link between Gaussian fields and Gaussian Markov random fields: The SPDE approach (with discussion). In: Journal of the Royal Statistical Society, Series B 73.4 (2011), pp. 423–498.
H. Bakka, H. Rue, G. A. Fuglstad, A. Riebler, D. Bolin, J. Illian, E. Krainski, D. Simpson, and F. Lindgren. Spatial modelling with R-INLA: A review. In: WIREs Computational Statistics 10:e1443.6 (2018). (Invited extended review). DOI: 10.1002/wics.1443.
E. T. Krainski, V. Gómez-Rubio, H. Bakka, A. Lenzi, D. Castro-Camilio, D. Simpson, F. Lindgren, and H. Rue. Advanced Spatial Modeling with Stochastic Partial Differential Equations using R and INLA. Github version . CRC press, Dec. 20

SPDE models

We call spatial Markov models defined on a mesh SPDE models.

SPDE models have 3 parts

A mesh
A range parameter $\kappa$
A precision parameter $\tau$

We use the SPDE effect to model the intensity of a point process that represents the locations of animal sightings.
Often such sightings are made by observers who cannot detect all the animals
To accurately estimate abundance, we require an estimate of the number of animals that remained undetected.

Thinned Point Process

The LGCP is a flexible approach that can include spatial covariates to model the mean intensity and a mean-zero spatially structured random effect to account for unexplained heterogeneity not captured by the covariates.

To account for the imperfect detection of points we specify a thinning probability function \[g(s) = \mathbb{P}(\text{a point at s is detected}|\text{a point is at s})\]
A key property of LGCP is that a realisation of a point process with intensity $\lambda(s)$ that is thinned by probability function $g(s)$, follows also a LGCP with intensity:

\[ \underbrace{\tilde{\lambda}(s)}_{\text{observed process}} = \underbrace{\lambda(s)}_{\text{true process}} \times \underbrace{g(s)}_{\text{thinning probability}} \]

Thinned Point Process

Lets visualize this on 1D: Intensity function with points

Thinned Point Process

Intensity (density) function with points and transect locations

Thinned Point Process

Detection function $\color{red}{g(s)}$
Here $\color{red}{g(s) =1}$ on the transects (at x = 10,30 and 50).

Thinned Point Process

Detection function $\color{red}{g(s)}$ and detected points

Thinned Point Process

The detection function describes the probability $\color{red}{p(s)}$ that an point is detected

Thinned Point Process

Observations are from a thinned Poisson process with intensity $\lambda(s) \color{red}{p(s)}$

Detection Function

Standard distance sampling approaches specify $g(s)$ as a function that declines with increasing distance
- horizontal distance to the observer for point transects
- perpendicular distance to the transect line for line transects

Detection Function

Standard distance sampling approaches specify $g(s)$ as a function that declines with increasing distance
- horizontal distance to the observer for point transects
- perpendicular distance to the transect line for line transects

Detection Function

Standard distance sampling approaches specify $g(s)$ as a function that declines with increasing distance
- horizontal distance to the observer for point transects
- perpendicular distance to the transect line for line transects
The thinning probability function is specified as a parametric family of functions.

Half-normal: $g(\mathbf{s}|\sigma) = \exp(-0.5 (d(\mathbf{s})/\sigma)^2)$

Hazard-rate :$g(\mathbf{s}|\sigma) = 1 - \exp(-(d(\mathbf{s})/\sigma)^{-1})$

Detection Function

Standard distance sampling approaches specify $g(s)$ as a function that declines with increasing distance
- horizontal distance to the observer for point transects
- perpendicular distance to the transect line for line transects
The thinning probability function is specified as a parametric family of functions.
The thinned-LGCP likelihood is given by:

\[ \pi(\mathbf{s_1},\ldots,\mathbf{s_m}) = \exp\left( |\Omega| - \int_{\mathbf{s}\in\Omega}\lambda(s)g(s)\text{d}s \right) \prod_{i=1}^m \lambda(\mathbf{s}_i)g(\mathbf{s}_i) \]

To make $g(s)$ and $\lambda(s)$ identifiable, we assume intensity is constant with respect to distance from the observer.
- In practice this means we assume animals are uniformly distributed with respect to distance from the line

Putting all the pieces together

The true point pattern $Y = \mathbf{s}_1,\ldots,\mathbf{s}_n$ are a realization of a Point process with intensity $\lambda(s)$

Putting all the pieces together

The true point pattern $Y = \mathbf{s}_1,\ldots,\mathbf{s}_n$ are a realization of a Point process with intensity $\lambda(s)$
We design a sampling survey to collect the data along transects

Putting all the pieces together

The true point pattern $Y = \mathbf{s}_1,\ldots,\mathbf{s}_n$ are a realization of a Point process with intensity $\lambda(s)$
We design a sampling survey to collect the data along transects
detected points are generated from the thinned PP with intensity $\color{red}{\tilde{\lambda}(s)}= \lambda(s)\color{red}{g(d(s))}$
- The log intensity $\log \color{red}{\tilde{\lambda}(s)} = \overbrace{\log \lambda (s)}^{\mathbf{x}'\beta + \xi(s)} + \overbrace{\log \color{red}{g(d(s))}}^{-0.5~d(\mathbf{s})^2\sigma^{-2}}$

Putting all the pieces together

The true point pattern $Y = \mathbf{s}_1,\ldots,\mathbf{s}_n$ are a realization of a Point process with intensity $\lambda(s)$
We design a sampling survey to collect the data along transects
detected points are generated from the thinned PP with intensity $\color{red}{\tilde{\lambda}(s)}= \lambda(s)\color{red}{g(d(s))}$
- The log intensity $\log \color{red}{\tilde{\lambda}(s)} = \overbrace{\log \lambda (s)}^{\mathbf{x}'\beta + \xi(s)} + \overbrace{\log \color{red}{g(d(s))}}^{-0.5~d(\mathbf{s})^2\sigma^{-2}}$
- The encounter rate, i.e. the number of observed animals within a distance $W$ follows $m \sim \text{Poisson} \left(\int_0^W \tilde{\lambda}(d)\text{d}d\right)$

The pdf of detected distances is $\pi(d_1,\ldots,d_m|m) \propto \prod_{i=1}^m\dfrac{\tilde{\lambda}(d_i)}{\int_0^W \tilde{\lambda}(d)\text{d}d}$

Putting all the pieces together

The true point pattern $Y = \mathbf{s}_1,\ldots,\mathbf{s}_n$ are a realization of a Point process with intensity $\lambda(s)$
We design a sampling survey to collect the data along transects
detected points are generated from the thinned PP with intensity $\color{red}{\tilde{\lambda}(s)}= \lambda(s)\color{red}{g(d(s))}$
- The log intensity $\log \color{red}{\tilde{\lambda}(s)} = \overbrace{\log \lambda (s)}^{\mathbf{x}'\beta + \xi(s)} + \overbrace{\log \color{red}{g(d(s))}}^{-0.5~d(\mathbf{s})^2\sigma^{-2}}$
- The encounter rate, i.e. the number of observed animals within a distance $W$ follows $m \sim \text{Poisson} \left(\int_0^W \tilde{\lambda}(d)\text{d}d\right)$

The pdf of detected distances is $\pi(d_1,\ldots,d_m|m) \propto \prod_{i=1}^m\dfrac{ g(d_i)}{\int_o^W g(d) \text{d}d}$ if $\color{red}{\tilde{\lambda}(d_i)} = \lambda \color{red}{g(d_i)}$

An approximation: Strips as lines

If the strips width ( $2W$ ) is narrow compared to study region ($\Omega$) we can treat them as lines.
- We need to adjust the intensity at a point $\mathbf{s}$ on the line to take account of the actual width of the strip
- Adjust the thinning probability to account for having collapsed all points onto the line.

An approximation: Strips as lines

The intensity at a point $\mathbf{s}$ on the line becomes $2W\lambda(s)$ instead of $\lambda(s)$.

Let $\pi(d)$ be the probability that the point is at a distance $d$ from the line.
Let $p(d)$ be the probability that is detected given it is at $d$.

Then, the thinning probability becomes $\pi(d)\times p(d)$, assuming the points are uniformly distributed within the strip then $\pi(d) = 1/W$ (the density of distances is assumed to be constant on the interval $[0,W]$).

This updates our thinning intensity to

\[ \log \tilde{\lambda}(s) = \underbrace{\mathbf{x}'\beta + \xi(s)}_{\log \lambda(s)} + \log p(d) + \log (1/W) \]

Typically $p(d)$ is a non-linear function, that is where inlabru can help via a Fixed point iteration scheme (further details available in this vignette)

Example

Example: Dolphins in the Gulf of Mexico

In the next example, we will explore data from a combination of several NOAA shipboard surveys conducted on pan-tropical spotted dolphins in the Gulf of Mexico.

A total of 47 observations of groups of dolphins were detected. The group size was recorded, as well as the Beaufort sea state at the time of the observation.
Transect width is 16 km, i.e. maximal detection distance 8 km (transect half-width 8 km).

Step 1: Define the SPDE representation: The mesh

First, we need to create the mesh used to approximate the random field. We can either:

Create a nonconvex extension of the points using the fm_mesh_2d and fm_nonconvex_hull functions from the fmesher package:

library(fmesher)

boundary0 = fm_nonconvex_hull(mexdolphin$points,convex = -0.1)

mesh_0 = fm_mesh_2d(boundary = boundary0,
                          max.edge = c(30, 150), # The largest allowed triangle edge length.
                          cutoff = 15,
                          crs = fm_crs(mexdolphin$points))

max.edge for maximum triangle edge lengths
cutoff to avoid overly small triangles in clustered areas

Step 1: Define the SPDE representation: The mesh

First, we need to create the mesh used to approximate the random field. We can either:

Use a pre-define sf boundary and specify this directly into the mesh construction via the fm_mesh_2d function

library(fmesher)


mesh_1 = fm_mesh_2d(boundary = mexdolphin$ppoly,
                    max.edge = c(30, 150),
                    cutoff = 15,
                    crs = fm_crs(mexdolphin$points))

max.edge for maximum triangle edge lengths
cutoff to avoid overly small triangles in clustered areas

Step 1: Define the SPDE representation: The mesh

All random field models need to be discretised for practical calculations.
The SPDE models were developed to provide a consistent model definition across a range of discretisations.
We use finite element methods with local, piecewise linear basis functions defined on a triangulation of a region of space containing the domain of interest.
Deviation from stationarity is generated near the boundary of the region.
The choice of region and choice of triangulation affects the numerical accuracy.

Step 1: Define the SPDE representation: The mesh

Too fine meshes $\rightarrow$ heavy computation
Too coarse mesh $\rightarrow$ not accurate enough

Step 1: Define the SPDE representation: The mesh

Some guidelines

Create triangulation meshes with fm_mesh_2d():
edge length should be around a third to a tenth of the spatial range
Move undesired boundary effects away from the domain of interest by extending to a smooth external boundary:
Use a coarser resolution in the extension to reduce computational cost (max.edge=c(inner, outer)), i.e., add extra, larger triangles around the border

Step 1: Define the SPDE representation: The mesh

Use a fine resolution (subject to available computational resources) for the domain of interest (inner correlation range) and avoid small edges ,i.e., filter out small input point clusters (0 $<$ cutoff $<$ inner)
Coastlines and similar can be added to the domain specification in fm_mesh_2d() through the boundary argument.
simplify the border

Step 1: Define the SPDE representation: The SPDE

We use the inla.spde2.pcmatern to define the SPDE model using PC priors through the following probability statements

$P(\rho < 50) = 0.1$
$P(\sigma > 2) = 0.1$

spde_model =  inla.spde2.pcmatern(
  mexdolphin$mesh,
  prior.sigma = c(2, 0.1),
  prior.range = c(50, 0.1)
)

Step 2: Define the Detection function

We start by plotting the distances and histogram of frequencies in distance intervals.

Then, we need to define a half-normal detection probability function. This must take distance as its first argument and the linear predictor of the sigma parameter as its second:

# define detection function
hn <- function(distance, sigma) {
  exp(-0.5 * (distance / sigma)^2)
}

Example: Dolphins in the Gulf of Mexico

The LGCP Model

\[ \begin{aligned} p(\mathbf{y} | \lambda) & \propto \exp \left( -\int_\Omega \lambda(\mathbf{s}) p(\mathbf{s}) \mathrm{d}\mathbf{s} \right) \prod_{i=1}^n \lambda(\mathbf{s}_i) p(\mathbf{s}_i)) \\ \eta(s) & = \color{#FF6B6B}{\boxed{\beta_0}} + \color{#FF6B6B}{\boxed{ \omega(s)}} + \color{#FF6B6B}{\boxed{ \log p(s)}} + \color{#FF6B6B}{\boxed{ \log (1/W)}}\\ \end{aligned} \]

The code

# define model component
cmp = ~ Intercept(1) + 
  space(main = geometry, model = spde_model) +
  sigma(1,
    prec.linear = 1,
    marginal = bru_mapper_marginal(qexp, pexp, dexp, rate = 1 / 8)
  ) 

# define model predictor
eta  = geometry + distance ~ space +
  log(hn(distance, sigma)) +
  Intercept + log(2)

# build the observation model
lik = bru_obs("cp",
              formula = eta,
              data = mexdolphin$points,
              ips = ips)

# fit the model
fit = bru(cmp, lik)

The integration scheme

# build integration scheme
distance_domain <-  fm_mesh_1d(seq(0, 8,
                              length.out = 30))
ips = fm_int(list(geometry = mesh,
                  distance = distance_domain),
             samplers = mexdolphin$samplers)

Example: Dolphins in the Gulf of Mexico

The LGCP Model

\[ \begin{aligned} p(\mathbf{y} | \lambda) & \propto \exp \left( -\int_\Omega \lambda(\mathbf{s}) p(\mathbf{s}) \mathrm{d}\mathbf{s} \right) \prod_{i=1}^n \lambda(\mathbf{s}_i) p(\mathbf{s}_i)) \\ \color{#FF6B6B}{\boxed{\eta(s)}} & = \color{#FF6B6B}{\boxed{\beta_0 + \omega(s) + \log p(s) + \log (1/W)}}\\ \end{aligned} \]

The code

# define model component
cmp = ~ Intercept(1) + 
  space(main = geometry, model = spde_model) +
  sigma(1,
    prec.linear = 1,
    marginal = bru_mapper_marginal(qexp, pexp, dexp, rate = 1 / 8)
  ) 

# define model predictor
eta  = geometry + distance ~ space +
  log(hn(distance, sigma)) +
  Intercept + log(2)

# build the observation model
lik = bru_obs("cp",
              formula = eta,
              data = mexdolphin$points,
              ips = ips)

# fit the model
fit = bru(cmp, lik)

The integration scheme

# build integration scheme
distance_domain <-  fm_mesh_1d(seq(0, 8,
                              length.out = 30))
ips = fm_int(list(geometry = mesh,
                  distance = distance_domain),
             samplers = mexdolphin$samplers)

Example: Dolphins in the Gulf of Mexico

The LGCP Model

\[ \begin{aligned} \color{#FF6B6B}{\boxed{p(\mathbf{y} | \lambda)}} & \propto \exp \left( -\int_\Omega \lambda(\mathbf{s}) p(\mathbf{s}) \mathrm{d}\mathbf{s} \right) \prod_{i=1}^n \lambda(\mathbf{s}_i) p(\mathbf{s}_i)) \\ \eta(s) & = \beta_0 + \omega(s) + \log p(s) + \log (1/W)\\ \end{aligned} \]

The code

# define model component
cmp = ~ Intercept(1) + 
  space(main = geometry, model = spde_model) +
  sigma(1,
    prec.linear = 1,
    marginal = bru_mapper_marginal(qexp, pexp, dexp, rate = 1 / 8)
  ) 

# define model predictor
eta  = geometry + distance ~ space +
  log(hn(distance, sigma)) +
  Intercept + log(2)

# build the observation model
lik = bru_obs("cp",
              formula = eta,
              data = mexdolphin$points,
              ips = ips)

# fit the model
fit = bru(cmp, lik)

The integration scheme

# build integration scheme
distance_domain <-  fm_mesh_1d(seq(0, 8,
                              length.out = 30))
ips = fm_int(list(geometry = mesh,
                  distance = distance_domain),
             samplers = mexdolphin$samplers)

Example: Dolphins in the Gulf of Mexico

The LGCP Model

The code

# define model component
cmp = ~ Intercept(1) + 
  space(main = geometry, model = spde_model) +
  sigma(1,
    prec.linear = 1,
    marginal = bru_mapper_marginal(qexp, pexp, dexp, rate = 1 / 8)
  ) 

# define model predictor
eta  = geometry + distance ~ space +
  log(hn(distance, sigma)) +
  Intercept + log(2)

# build the observation model
lik = bru_obs("cp",
              formula = eta,
              data = mexdolphin$points,
              ips = ips)

# fit the model
fit = bru(cmp, lik)

The integration scheme

# build integration scheme
distance_domain <-  fm_mesh_1d(seq(0, 8,
                              length.out = 30))
ips = fm_int(list(geometry = mesh,
                  distance = distance_domain),
             samplers = mexdolphin$samplers)

Results

Results: posterior summaries

We can use the fit$summary.fixed and summary.hyperpar to obtain psoterior summarie sof the mdoel parameters.

	mean	0.025quant	0.975quant
Intercept	−8.41	−9.47	−7.62
sigma	−0.05	−0.46	0.36
Range for space	131.74	41.79	320.28
Stdev for space	1.17	0.72	1.78

The spde.posterior allow us to plot the posterior density of the Matern field parameters

spde.posterior(fit, "space", what = "range") %>% plot()

Results: posterior summaries

We can use the fit$summary.fixed and summary.hyperpar to obtain posterior summaries of the model parameters.

	mean	0.025quant	0.975quant
Intercept	−8.41	−9.47	−7.62
sigma	−0.05	−0.46	0.36
Range for space	131.74	41.79	320.28
Stdev for space	1.17	0.72	1.78

The spde.posterior allow us to plot the posterior density of the Matern field parameters

spde.posterior(fit, "space", what = "log.variance") %>% plot()

Results: predicted densities

To map the spatial intensity we first need to define a grid of points where we want to predict.

We do this using the function fm_pixel() which creates a regular grid of points covering the mesh
Then, we use the predict function which takes as input
- the fitted model (fit)
- the prediction points (pxl)
- the model components we want to predict (e.g., $e^{\beta_0 + \xi(s)}$)
To plot this you can use ggplot and add a gg() layer with your output of interest (E.g., pr.int$spatial)

library(patchwork)
pxl <- fm_pixels(mesh, dims = c(200, 100), mask = mexdolphin$ppoly)
pr.int <- predict(fit, pxl, ~ data.frame(spatial = space,
                                      loglambda = Intercept + space,
                                      lambda = exp(Intercept + space)))

ggplot() +
  gg(pr.int$spatial, geom = "tile")

Results: predicted densities

We can also use the predict function to predict the detection function:

distdf <- data.frame(distance = seq(0, 8, length.out = 100))
dfun <- predict(fit, distdf, ~ hn(distance, sigma))
plot(dfun)

Results: predicted expected counts

We can look at the posterior for the mean expected number of dolphins:

predpts <- fm_int(mexdolphin$mesh, mexdolphin$ppoly)
Lambda <- predict(fit, predpts, ~ sum(weight * exp(space + Intercept)))
Lambda

      mean       sd   q0.025     q0.5   q0.975   median sd.mc_std_err
1 320.8391 105.9736 193.2869 303.2979 614.2154 303.2979      10.92919
  mean.mc_std_err
1         12.7832

Results: predicted expected counts

We can also get Monte Carlo samples for the expected number of dolphins as follows:

Ns <- seq(50, 450, by = 1)

Nest <- predict(fit, predpts,
  ~ data.frame(
    N = Ns,
    density = dpois(
      Ns,
      lambda = sum(weight * exp(space + Intercept))
    )
  ),
  n.samples = 2000
)

Nest <- dplyr::bind_rows(
  cbind(Nest, Method = "Posterior"),
  data.frame(
    N = Nest$N,
    mean = dpois(Nest$N, lambda = Lambda$mean),
    mean.mc_std_err = 0,
    Method = "Plugin"
  )
)

Model comparison

Fitting a Model with a hazard detection function

hr <- function(distance, sigma) {
  1 - exp(-(distance / sigma)^-1)
}
eta_2 <- geometry + distance ~ space +
  log(hr(distance, sigma)) +
  Intercept + log(2)

lik_2 = bru_obs("cp",
              formula = eta_2,
              data = mexdolphin$points,
              ips = ips)

# fit the model
fit_hazard = bru(cmp, lik_2)

pr.int1 <- predict(fit_hazard, pxl, ~ data.frame(spatial = space,
                                      loglambda = Intercept + space,
                                      lambda = exp(Intercept + space)))
distdf <- data.frame(distance = seq(0, 8, length.out = 100))
dfun1 <- predict(fit_hazard, distdf, ~ hr(distance, sigma))

ggplot() +
  gg(pr.int1$loglambda, geom = "tile") +
  scale_fill_scico(palette="imola",name=expression(log(lambda)))+
  plot(dfun1)+plot_layout(ncol=1)

Model comparison with GoF

Look at the goodness-of-fit of the two models in the distance dimension

bc <- bincount(
  result = fit,
  observations = mexdolphin$points$distance,
  breaks = seq(0, max(mexdolphin$points$distance), length.out = 9),
  predictor = distance ~ hn(distance, sigma)
)
attributes(bc)$ggp

bc1 <- bincount(
  result = fit_hazard,
  observations = mexdolphin$points$distance,
  breaks = seq(0, max(mexdolphin$points$distance), length.out = 9),
  predictor = distance ~ hn(distance, sigma)
)
attributes(bc1)$ggp

Lecture

Outline

What is INLA? What is inlabru?

Where?

So… Why should you use inlabru?

Distance Sampling

Overview of Distance Sampling

Density surface models

Density surface models

Density surface models

Point Processes

Point process data

Defining a Point Process

Homogeneous Poisson Process

Inhomogeneous Poisson process

Inhomogeneous Poisson process

Inhomogeneous Poisson process

Limitations with IPP

The Log-Gaussian Cox Process

Gaussian Random Fields

Gaussian Random Fields

The SPDE approach

The SPDE approach

Penalized Complexity (PC) priors

Learning about the SPDE approach

SPDE models

Thinned Point Process

Thinned Point Process

Thinned Point Process

Thinned Point Process

Thinned Point Process

Thinned Point Process

Thinned Point Process

Thinned Point Process

Thinned Point Process

Thinned Point Process

Detection Function

Detection Function

Detection Function

Detection Function

Putting all the pieces together

Putting all the pieces together

Putting all the pieces together

Putting all the pieces together

Putting all the pieces together

An approximation: Strips as lines

An approximation: Strips as lines

Example

Example: Dolphins in the Gulf of Mexico

Step 1: Define the SPDE representation: The mesh

Step 1: Define the SPDE representation: The mesh

Step 1: Define the SPDE representation: The mesh

Step 1: Define the SPDE representation: The mesh

Step 1: Define the SPDE representation: The mesh

Step 1: Define the SPDE representation: The mesh

Step 1: Define the SPDE representation: The SPDE

Step 2: Define the Detection function

Example: Dolphins in the Gulf of Mexico

Example: Dolphins in the Gulf of Mexico

Example: Dolphins in the Gulf of Mexico

Example: Dolphins in the Gulf of Mexico

Results

Results: posterior summaries

Results: posterior summaries

Results: predicted densities

Results: predicted densities

Results: predicted expected counts

Results: predicted expected counts

Model comparison

Fitting a Model with a hazard detection function

Model comparison with GoF

What is INLA? What is `inlabru`?

So… Why should you use `inlabru`?